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The transport and fluctuation properties of organic molecules ordered parallel between two Au 
contact leads are investigated by the method of surface Green function. From first-principles sim- 
ulation the relevant hopping parameters are extracted and used to calculate nonlinear transport 
coefficients with respect to an external bias voltage. A staggering of conductance is found in de- 
pendence on the number of molecules squeezed in-between the contacts. The thermal properties 
show an anomalous behavior whenever the voltage reaches the values of the molecular energy levels 
active for transport. The thermoelectric figure of merit shows a resonance allowing to reach values 
even larger than one. 
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I. INTRODUCTION 

The goal to develop low-cost storage and microelec- 
tronic devices has triggered an enormous activity in the 
research of organic field effect transistors (OFETs) based 
on different polymers 1 -!^ and small organic molecules^. 
It is desirable that the molecular material possesses a 
high structural ordering 5 ' 6 in order to reach high charge 
carrier mobilities and low resistive losses. Among the 
most promising materials are oligothiophenes and their 
derivatives 7 -. This is due to the variety of intra- and in- 
termolecular interactions, which originates from the po- 
larizability of the sulfur electrons and the embedding aro- 
matic 7r-electron system&Sii .. 




FIG. 1: The thiophene molecules between 2 Au contacts. 
The performance of OFETs based on oligothiophenes 3 



shows characteristic features, e.g. the current starts at 
a certain threshold of gate voltage and reaches a satu- 
ration value for certain drain voltages. In general, the 
charge transport occurs in the direction perpendicular to 
the plane of the thiophene rings due to good 7r-stacking. 
This suggests to construct devices with parallelly or- 
dered molecular rings. This is the perpendicular direc- 
tion to the usually considered transport through vari- 
ous molecule classes ranging from transport dependent 
on the thickness 1 ^ to metallic behaviour independent of 
the thickness of molecules^. Thick-film devices based on 
thiophenes were reported several years ago^ 3 -. Recently 
an OFET structure has been built from ultra-thin self- 
assembled films made up from oligothiophenes, which are 
arranged in a highly-order lamellar stacking perpendic- 
ular to the substrate surface 3 - as illustrated in figures Q] 
and [3] 

The energy gap between the lowest unoccupied 
(LUMO) and the highest occupied molecular level 
(HOMO) measures about 3 eV in quarterthiophene. 
For the molecular structure see figure [TJ In contrast 
to the conduction and valence band in semiconductors 
the transport is due to localized states dominated by 
hoppin g 14 i 15 rather than due to delocalized states lim- 
ited by the scattering in semiconductors. 

If one considers the charge transport between the 
molecules one can think of two extreme mechanisms. On 
one side one has either the ballistic transport through 
bonds, or the hopping transport through space. The dif- 
fusive transport on the other side is present if the car- 
riers meet scattering partners. For the present situation 
of parallely stacked molecules the hopping transport by 
tunneling between the localized HOMO or LUMO states 
is the relevant one. A related mechanism of charge trans- 
fer due to the bending motion of molecules by touching 
might lead to shuttling transport 16 which could be in- 
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FIG. 2: The parallel stacking of thiophene molecules be- 
tween 2 Au contacts building the source (S) and drain (D) 
below. Typical sizes are W = 1000/xm, T — 20nm and 
L = 10/im. One slice is enlarged (above) presenting the chem- 
ical structure-. 



vestigated by resonant spectroscopy^. Of course, the 
most crucial question is the tilting of molecules triggered 
by the metal-molecul contacts^ which requires a refined 
interface engineering-^. 

Here we construct a hopping transport picture and 
want to investigate the particle and thermal transport 
properties as well as current fluctuations. We will project 
the electron transport including Coulomb and inelastic 
tunneling to a simple one-particle hopping Hamiltonian 
by calculating the matrix elements through overlap inte- 
grals of orbitals obtained from density functional calcu- 
lation. We account for the Coulomb force only in deter- 
mining these effective hopping elements. The Coulomb 
interaction is not considered explicitly in the hopping 
Hamiltonian such that we concentrate only on coherent 
motion but under nonlinear external bias voltage. The 
Coulomb effect to the hopping transport has been dis- 
cussed e.g. in [20.21.22]. One might expect that the 
Coulomb effect will lead to a blocking of the transport 
due to the creation of a Coulomb gap. This is beyond the 
scope of the paper, and since we deal with a low charge 
per molecule of ~ 10 -6 , we neglect this effect at room 
temperature. Other inelastic tunneling processes we do 
not account explicitly in the hopping model than in the 
cumulative sense of DFT based hopping matrix elements. 



Therefore this model is an effective one and cannot claim 
to resolve detailed information on the single molecular 
level. The predictive power of such an effective model 
can finally be checked only with corresponding experi- 
mental data which are expected in nearest future. 

In calculating the transport properties for hop- 
ping situation the Greenfunctions are a convenient 
too l 2 ?i 24 ' 2 ?i 2 9i 2 Ti 2 9 . We will use the method of surface 
Greenfunctions 29,30 . Alternative methods of direct nu- 
merical inversions or expansions with the help of the Flo- 
quet theorem can be found in the literature, e.g£^. The 
method of surface Greenfunctions has the advantage to 
be numerically fast (linear proportional to site numbers) 
and free of large compensations. Furthermore it allows to 
describe a time-dependent response in a convenient way 29 
and can be easily extended to multiband problems 3 - 2 -. 

First we give the basic ideas of the surface Greenfunc- 
tion method in chapter II with the approach to the cur- 
rent and current fluctuations. The technical details of the 
numerically extremely fast method are presented in the 
appendix. In chapter III we describe the first-principles 
approach to derive the necessary hopping parameters. 
The numerical results are presented in chapter IV in- 
cluding the nonlinear transport coefficients leading to an 
anomaly in the thermal transport. The summary dis- 
cusses possible experimental verifications and the status 
of current theoretical developments. 

II. METHOD OF SURFACE GREENFUNCTION 

A. Convenient formulation of the boundary 
problem 

We consider the tight-binding Hamiltonian 

H= £ (\j)v j ij\ + \j+l)t j+l ij\ + \j-l)t j ij\\ (1) 

3=—ao V 

describing N molecules in-between the crystal leads at 
the left, j < 0, and the right side, j > N + 1. The energy 
levels are Vj and the hopping between neighbouring sites 
are tj. 

For the formulation of the scattering one might be 
tempted to use the Lippmann-Schwinger formulation. To 
this end one writes the Hamiltonian as a sum of the un- 
perturbed crystal, H°, made by an extension of the left 
semiinfinte crystal all over the whole space, and the per- 
turbation, H' = H — H . An incoming wave ip is not 
yet influenced by the barrier of the leads, thus it is the 
eigenstate of H°, 

H°tP° = EiP°. (2) 

Eigenstates of the homogeneous system are plane waves 
with a momentum k related to the energy by E = 
2vq — 2vq cos(ka) with the distance of atoms a. The 
wave approaching the barrier is ip°(x) = e-xjp(ikx) or 
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discretized xj — j a in the tight-binding representation 
ipj = exp(ikaj). 

The total wave function -0° + ip' solves Schrddinger's 
equation 

(H° + H')(iP°+iP') = E(ip +iP'), (3) 

which can be rearranged with the help of ([2]) into the 
Lippmann-Schwinger equation, 

(E - H)ip' = H'ip°. (4) 

The boundary condition that the diffracted wave prop- 
agates from the barrier to infinity is realized by an in- 
finitesimal shift of the energy E into the complex plane. 

The above straightforward implementation of the 
Lippmann-Schwinger idea becomes numerically ex- 
tremely inconvenient if the left and right leads are dif- 
ferent by hopping, energy levels and external bias. For 
instance, when the tunneling junction is biased by the 
voltage V, the perturbation H' in the right half-space 
equals the potential of the bias voltage, H'-- — v N+ i for 
all j > N+l. The perturbation then extends everywhere 
except for the left lead, and the source term in the right 
hand side of (0J is nonzero over this region, too. 

In this case, the Lippmann-Schwinger separation of the 
wave function into an infinite plane wave ip° and correc- 
tion ip is not a favourable starting step. This is clearly 
seen for the case when electrons have an energy at the 
level of the right lead. All electrons are then reflected at 
the barrier and only exponential tails of the wave func- 
tion penetrate into the lead. The finite value of ip° has to 
be compensated by ip. This requires a large source term, 
H'lp , and a very accurate treatment of the Lippmann- 
Schwinger equation. In fact, only identities can guaran- 
tee a correct compensation of the large parts of incoming 
and outgoing waves. 

There is a simple modification of the Lippmann- 
Schwinger idea by which we can circumvent the pene- 
tration of the incoming wave into the barrier and the 
right lead. Let us cut the incoming wave by putting 

tpf = ip] for j < 0, and 

ip? = for j > 1. (5) 

The wave function is now split differently than before, 
tp° + ip 1 = i/j° + ip, and Schrddinger's equation, (E — 
H)(ip c + ip) = 0, gives 

{E-H)ip = Q, n = -(E-H)?p c . (6) 

The term Q is a well behaved source which is nonzero 
only at two layers, j = and j = 1, with 

= -ti$ - -he ika , Ox = tx^g - h, (7) 

which can be easily proven calculating (j \ E — H \ tp c ) 
with the help of ([1]). The wave function ip does not in- 
clude the incoming part, we can thus use the retarded 



boundary condition, z = E + iO, and write the wave 
function in terms of Greenfunction, 

%p = —^—tt = Gil (8) 

z — H 

or explicitly 

=G J -ifi 1 + Gj fi - (9) 



B. Reflection and transmission 

When considering the left lead, the incoming wave 
shows the reflected part characterized by the reflection 
coefficient r 

0_i = re ika 

= t SLi(Goifii + Goono) (10) 

where we used the recursion formula (IA22p for the Green- 
function. Since the energy conservation reads cos (fca) = 
(z — vo)/2to we employ the explicit from of the surface 
Greenfunction in the lead (|A15[) to derive the relation 

t S l = e* ka . (11) 

This allows us to identify from (|10p the reflection coeffi- 
cient 

r = G if!i+Goono. (12) 

Considering the right lead we have the transmission co- 
efficient according to 

xP N+1 = tc^+V 

— Gjv+i,i"i + Gn+i,oHo 

= tN +1 S 1 N+1 (GN t iili + Gn.o^o) (13) 

where we used the recursion relation (|A21|) . Since the 
energy conservation on the right lead means cos (n)a — 
{z — VN+i)/2tN + i we employ the explicit form of the 
surface Greenfunction in the lead (IA24p to obtain the 
relation 

tw + iS r N+1 = e* Ka (14) 
and from Q13[) follows the transmission coefficient 

t = e~ lKaN (Gn,i^i + Gno^o) 
= -2ihe- iKaN G Nfi sin (ka) 
= -2^ie- lKa ( Ar+1 )G w+1 , sin(fca) (15) 

representing different equivalent forms. The reflection 
R = \r\ 2 and transmission T = \t\ 2 obey T + R = 1, of 
course. 
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C. Current and current fluctuations 

For the tight-binding Hamiltonian (Q} the velocity op- 
erator x = j- [H, x] is easily computed from the position 

operator x = aJ^J \ j) 01 as 

i 



which is the known Landauer-Biittiker form of 
conductanc e 27 ' 31 . 

The time-averaged current autocorrelation function 3 ^ 



J 2 = lim - 



dt 



dr(AJ(t)AJ(t - t)) (22) 



i = f E^+i(i^o- + ii-b + i>oi)- as) 

3 

The mean current j = (ex) is conserved. Therefore it 
is sufficient to calculate the current which runs from the 
molecules to the lead on the right side 



lea 



(17) 



where we introduced the definition of the correlated 
Greenfunctions G < . The Greenfunction provided in 
the appendix are causal functions. According to the 
Langreth/Wilkins rules^ 3 . which establish the general- 
ized Kadanoff and Baym formalis m 34 ' 35 from the prod- 
uct of causal functions A and B one obtains the cor- 
related part via the rule [AB)< = A R B< + A<B A . 
The retarded and advanced functions obey the relation 
i(A R - A A ) = (A < - A<) for fcrmions. With the help 
of these rules and the recursion relation of the causal 
Greenfunction (|A21|) we can write 

= tN+i{S§ +1 G^ N + Sm +1 G nn ) (18) 



and repeatedly applying (|A13|) the current ([TTJ) reads 



N+l N-l 



3=1 



k=l 



Further simplification can be achieved by rewriting the 
diagonal Greenfunction Gnn in (|19p into the Greenfunc- 
tion Gjv+i,o with the help of (|A22|) . We observe that the 
surface correlated Greenfunction is the one of the left 
lead. It can be expressed in terms of the spectral func- 
tion Ai = i(S l R — S l Q A ) and the Fermi-Dirac function /; 
describing the occupation of the left lead via Sq = A l fi 
and Sq = A l (l — /;) and analogously for the right lead. 
One obtains finally for the current J = j /a 

J(z) = ^\G N+1 , (z)\Hlt% +1 A r N+1 (z)A l (z) 



fr 



z-eV 
T 



-fl 



In order to obtain the conductance we have to integrate 
over all energies and to divide by the voltage V = U/e 



Hik^^-'A < 21 > 



describes the zero-frequency noise of the current (|20p . 
The result reads^ 7 - 3 ! 



J 



2ttHU 



7jv+i,o(/r — fi) 



+ (1 - /r)/r + (1 - fl)fl 



(23) 



As a measure for the fluctuations we will use the Fano 
factor 



F 



eJ' 



(24) 



III. 



CALCULATION OF THE TIGHT-BINDING 
PARAMETERS 



Projecting the transport problem to a nearest neig- 
bour hopping we do need to calculate appropriate 



hopping parameter tj 



1 3,3+1- 



The electronic 



coupling matrix elements Tij are computed using a 
fragment orbital approach within the self-consistent 
charge density-functional based tight-binding (SCC- 
DFTB) method 3 ^ 3 ^^^! 4 ^ A similar approach has 
been already successfully applied to calculate charge- 
transfer matrix elements for hole transfer in DNA— 

The electronic coupling between two molecular orbitals 
<j)i and <j)j of neighboring molecules or molecular frag- 
ments i and j is calculated from isolated molecules by 
the hopping matrix elements 



Tj 



(25) 



where Hks is the Kohn-Sham Hamiltonian^. By apply- 
ing the atomic-basis-set expansion 



(26) 



the coupling integrals in the molecular-orbital basis can 
be efficiently evaluated as 



Cu4 (r} ll \H K s\'n v 



^c^H^. (27) 



(20) H^ v is the Hamilton matrix represented in the atomic 



basis. 



This matrix and the atomic coefficients c' are 

a* 



estimated using the SCC-DFTB method 3 ^^^ 4 ^ 3 .. 

The SCC-DFTB method is based on the density- 
functional theory of Hohenberg and Kohn 4 ! in the for- 
mulation of Kohn and Sham 4 ^. The single-particle Kohn- 
Sham eigenfunctions (pi are expanded in a set of localized 
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FIG. 3: The tight binding parameter for nearest neighbor 
hopping of thiophene versus the distance. The equilibrium 
distance is marked by a vertical line. The variance of the 
distance of molecules due to thermal motion is indicated as a 
box. 



atom-centered basis functions tj^ (Eq |2l))) . These func- 
tions are determined by self-consistent density-functional 
calculations on the isolated atoms employing a large set 
of Slater-type basis functions which obey 



■VeffiPo 



(28) 



The calculation of the Ty according to (|2"5|) is based on 
orthogonal basis functions, 



Spy = (r}fi\v^) = <W, 



(29) 



where the atomic basis functions at different atomic cen- 
ters have been orthogonalized using the Schmidt orthog- 
onalization 



to?) 



Hence, the Hamilton matrix 



VlJ.WKS-DFTB\fjv 



(30) 



(31) 



can be calculated for all necessary combinations of or- 
bitals f]^ and f] v on atoms a and /3, and stored in tables. 

First, the isolated molecule i has been geometry- 
optimized and the atomic coefficients c* of the respec- 
tive orbital fj, have been estimated by standard DFTB 
calculation o 38 i 39 ' 40 i 41 ' 42 i 43 . Second, the atomic Hamilton 
matrix (f3"Tj) is constructed in order to calculate the lo- 
using the orthogonalized and non-confined atomic basis 
functions rj^. 

The Tij have been computed between the highest oc- 
cupied orbitals (HOMOs) and lowest unoccupied orbitals 
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FIG. 4: The work function of thiophene molecules with re- 
spect to Au contacts. 



(LUMOs) of two quarterthiophene molecules. The re- 
spective nearest neighbour hopping parameters tj — 



T 



3,3+1 



in dependence on the stacking distance of the 



molecules is shown in figure [3l 



IV. RESULTS 
A. Conductance for thiophene molecules 

We present now the results for thiophene molecules 
between Au contacts. The work function of Au is —5.1 
eV and the HOMO levels of thiophene are at —5.5 eV 
such that Vj = —0.4 eV as illustrated in figure |4j The 
LUMO levels are much higher at —3.8 eV and do not 
play any role for the conductance here. The calculated 
hopping parameters of the last section vary slightly with 
the distance of the molecules as indicated by the box in 
figure [3] This translates into a slight variation of the 
conductance as well. We have plotted in figure [5] the 
dependence of the linear conductance at infinitesimal bias 
on the hopping parameter. A variance of the distance 
of the molecules between 2 and 4 A results only in a 
variation of the conductance by around 10% due to the 
flat curves in figures [3] and 

It is interesting to discuss the dependence of the 
conductance on the number of molecules squeezed in- 
between the contacts. In figure[|5]we plot the conductance 
as a function of the number of molecules for different ap- 
plied voltages. One recognizes a staggering up to 20-25 
molecules. Above this critical number the conductance 
saturates. With increasing applied voltage this stagger- 
ing is already damped out for a smaller number of par- 
ticles besides the overall decreasing of the conductance 
which will be discussed below. 

In order to find the origin of this effect let us plot the 
conductance for different tight binding parameter. In fig- 
ure [7] we see that the overall conductance becomes larger 
for larger tight binding parameter which is obvious due 
to the better coupling between the molecules. Moreover 
one can see that the period of staggering becomes shorter 
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FIG. 5: The conductance of 20 thiophene molecules between 
Au contacts versus the tight binding parameter. The variance 
according to figure [3] is indicated by the box. 



FIG. 7: The conductance versus the number of thiophene 
molecules between Au contacts for different tight binding pa- 
rameter in the linear response regime of applied voltage. 
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FIG. 6: The conductance versus the number of thiophene 
molecules between Au contacts for different applied voltages. 



for larger tight binding parameter. This suggests that we 
possibly see here a coherence effect over a certain number 
of molecules which is controlled by the hopping parame- 
ter. 



B. Density of states, current and current 
fluctuation 

The conductance is sensitively dependent on the rela- 
tive hopping parameter in the Au lead compared to the 
molecule ones. In the next figure [8] the hopping parame- 
ter in the Au leads are varied and the nonlinear conduc- 
tance is plotted versus applied voltage bias. 



FIG. 8: The conductance versus applied voltage bias for dif- 
ferent hopping parameters from the lead to the thiophene. 
The energy levels of the leads are at vo = and the HOMO 
energy levels of the thiophene are at Vj = —0.4. As refer- 
ence curve the case of homogeneous lead material, Vj = 0, is 
plotted as well (thin line). 



Compared with the homogeneous case of the lead ma- 
terial, vq = Vj = 0, which gives the maximal conduc- 
tance, one sees that the mismatch between the lead hop- 
ping parameter and the molecule hopping parameter low- 
ers the conductance. While a smaller hopping parameter 
in the leads lead to an overall shrinkage of the current- 
voltage curve, a large hopping parameter in the lead leads 
to a lower but broader curve. This means that the de- 
vice is conductive even for higher applied voltages where 
the homogeneous case does not allow a current any more 
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FIG. 9: The density of states according to the cases in figure[8] 
The allowed region for transmission is indicated as thick line 
on the x-axes. 



since the energy dispersions of the left and right layer do 
not overlap. This is of course a mechanism which can be 
used as molecular switc h 48 i 49 . 

The different conductance can be understood from the 
density of states which is simply given by the diagonal 
Grccnfunction 



N+l 

D N = ImGj 

3=0 



(32) 



In figure [9] we give the density of states according to fig- 
ure [5J One recognizes the appearance of bound and res- 
onance states below the conductance range for lead hop- 
ping parameters less than or equal to the molecular ones, 
(b) and (c). For larger hopping parameters, (d)-(f), the 
bound states disappear. This can be easily understood 
since electrons are more easily hopping from the leads 
to the molecules than transported between the molecules 
which prevents any standing wave or bound state. 

Next we calculate the corresponding current-current 
fluctuations via the Fano factor (|2"4")1 . In figure [TO] we 
plot the Fano factors corresponding to the situation of 
figure [HI One recognizes that the fluctuations are sup- 
pressed with increasing lead to molecule hopping. For 
larger bias voltage, they approach the Schott noise limit, 
i.e. Poissonian characterized by F = 1. It is remark- 
able that the noise at zero bias becomes super- Poissonian 
for hopping parameters of the lead smaller than the 
molecules. For higher hopping parameter these fluctu- 
ations are suppressed again. 

As a next question we consider the optimal hopping 
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FIG. 10: The Fanofactor versus applied voltage bias for dif- 
ferent hopping parameters as used in figure [8] 
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FIG. 11: The conductance G versus applied voltage bias for 
different barrier height of the leads. The hopping parameter 
is to = 0.4 and the thiophene HOMO levels relativ to the 
leadswj = —0.4. 



parameter between Au contacts and the molecules but 
vary the barrier of the leads. This is motivated by the 
barrier that occurs at charge injection at metal-organic 
interfaces described by the image potential 50 . In figurefTTI 
the conductance is plotted for different forms of the lead 
barrier. 

With increasing barrier height the conductance is low- 
ered as one would expect. The corresponding fluctuations 
in figure [TO] show that the fluctuations are suppressed 
with decreasing barrier. This means that the quality of 
the transport can be improved if the lead to molecule 
barrier is fabricated as low as possible. 

The corresponding density of states in figure [TO] show 
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C. Transport coefficients 
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With the same tools as used so far it is very easy to 
compute all transport coefficients, even bias voltage de- 
pendent. We assume a different temperature on the right 
lead T r = T + AT/2 and on the left lead T r = T — AT/2. 
Linearizing with respect to the temperature gradient, but 
not with respect to the voltage, the particle current (|20|) 
reads 



J = L n V + L 12 AT 
with the nonlinear Onsager coefficients 

'z-eV 



(33) 



L 11 = 



-1 1 
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FIG. 12: The Fanofactor versus applied voltage bias for dif- 
ferent barrier heights of the leads as used in figure 1111 
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fr rp2 fl 



(34) 



with /' = — /(l — /). With absent temperature gradients 
we have the conductance (JUJ) 



G = L 



li 



(35) 




tr z feV] 



The voltage compensating the current due to the tem- 
perature gradient, V = —QAT, such that J = deter- 
mines the thermoelectric power (Seebeck coefficient) as 



(a) vo = 0.4 



.J -J.S 




-TT„z [eV] 



Q 



L 11 ' 



(36) 
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(d) v = -0.2 



Analogously we can give the heat current J q which is 
just (|20[) but with an additional energy factor z under the 
integrand. With small temperature gradients the heat 
current takes the form 



(b) vo = 0.2 
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J q = L V + AT 
with the nonlinear Onsager coefficients 

- eV 



(37) 



(e) Vo = -0.4 



L 21 

L 22 



e 
W 



e 

'2h 



dzzT N+1 . 
dzzT N+1:0 



fr 



T 

eV 



y2 



fr rp2 fl 



(38) 



(c) no = 

FIG. 13: The density of states according to the figures in 
figure [TT1 The allowed region for transmission is indicated as 
thick line on the x-axes. 



that the highest lead barrier (a) yields the largest num- 
ber of bound states below the conductance range. When 
lowering the barrier the bound state disappears and the 
density of states for the almost homogeneous case (e) 
appears. The remaining difference to the homogeneous 
case comes from the level mismatch of the HOMO and 
the leads. 



The thermoelectric conductance is measured if we keep 
the particle current zero, J = 0, which results from 
and (|37p into 



J q = KAT 
with the thermoelectric conductance 

L 21 L 12 



K = L 22 - 



L 



li 



(39) 



(40) 



The ratio between the thermoelectric conductance and 
the conductance is called Wiedemann-Franz law and 
should be linearly proportional to the temperature 



K 
G 



^22 



L 12 L 2 



(41) 
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FIG. 14: The conductance G, thermal heat conductance K, 
Wiedemann Franz law, thermopower Q and figure of merit 
versus temperature for different applied voltage bias near the 
linear response. 



FIG. 15: The conductance, heat conductance, Wiedemann 
Franz law, thermopower and figure of merit versus temper- 
ature for different applied voltage bias far beyond linear re- 
sponse. 



As a measure for the effectiveness of the thermoelectric 
devices the dimensionless figure of merit or ZT factor is 
often presented^ as 

ZT T Q 2 G T(L 12 ) 2 

K L n L 22 - L 21 L 12 V ' 

which has become an important quantity for thin film 
thermoelectric device o 52 ' 53 . 

In figures [T4land[T5lwe present the results for the trans- 
port coefficients dependent on temperature and the ex- 
ternal bias. In figure [14] we plot the results for low ex- 
ternal bias to show at which voltage one sees deviations 
from the linear response result. With increasing applied 
voltage the conductance is lowered since the overlap be- 
tween left and right conductance channel shrinks. This 
was seen already from figures [51 andfTTl The thermal con- 
ductance shows a nonlinear temperature behavior devel- 
oping a maximum at certain temperatures. This shows 
that the thermal conduction is much more sensitive to 
the interplay between external voltage and HOMO lev- 



els of the molecules. The thermopower reveals a similar 
sensitivity. The Wiedemann Franz law is seen to be not 
fullfillcd strictly. Instead we observe a deviation from the 
linear temperature behaviour of up to 10% for lower tem- 
peratures. The figure of merit for low applied voltages is 
negligible but plotted for completeness. 

The situation is changing if we apply higher voltages 
as shown in figure 1151 The observed maximum in the 
thermal conductance becomes more pronounced and for 
low temperatures the thermal conductance even changes 
sign if we approach an external voltage comparable to 
the HOMO levels. This has dramatic consequences on 
the figure of merit. At the temperature where the ther- 
mal conductance is changing sign we observe a resonance 
structure enhancing the figure of merit dramatically. For 
applied voltages equal to the HOMO levels we have the 
largest resonance but at temperatures around 80K. The 
resonance is shifted towards higher temperatures for volt- 
ages slightly above and below, but becomes broader, as 
well. Interestingly it can be observed that the ther- 
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mopower is changing sign for certain temperatures and 
becomes completely negative for voltages equal to the 
HOMO levels. This shows that for such voltages the 
thermal current is reversing sign which can be seen anal- 
ogously to the refrigerator effect described in [54| by an 
energy-selective transmission of electrons through a spa- 
tially asymmetric resonant structure subject to ac driv- 
ing. 



V. SUMMARY 

Organic molecules (thiophene) ordered paralel between 
Au contacts are considered under external voltage bias. 
The charge transport due to hopping is investigated and 
the particle and thermal transport properties calculated 
with the help of the numerically very fast surface Green- 
function technique. A tight binding Hamiltonian is used 
as underlying model where the hopping parameter are 
extracted from density functional calculations. 

The conductance as well as the current fluctuations 
can be shaped by changing the lead barrier height. This 
is due to the bound and resonance states which can be 
created with high enough barriers. The nonlinear con- 
ductance due to high voltage bias is discussed as well 
in dependence on different hopping parameter which de- 
pend strongly on the stacking distance of the molecules. 

We report a staggering of the conductance dependent 
on the number of molecules squeezed inbetween the two 
contacts. This finite-number-of-molecule effect we traced 
down to a coherence phenomenon controlled by the hop- 
ping parameter. 

The thermal conductance and the thermopower shows 
a nonlinear behaviour with higher applied voltages even 
changing sign at special temperatures. This results into a 
resonance structure in the figure of merit and values near 
unity. The latter effect was found to be dependent on the 
applied voltage most pronounced if the voltages is of the 
order of the energy levels of the HOMO relative to the 
leads. Besides application of these organic molecules in 
field effect transitor like structures, these materials might 
be of interest also for thermoelectric elements. 
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APPENDIX A: MATRIX INVERSION WITH 
SURFACE GREENFUNCTIONS 

The method of surface Grecnfunctions provides an ex- 
tremely fast numerical inversion method to solve the 
Greenfunction equation 



G(z -H) = l 



(Al) 



which means the inversion of the matrix z — H. Here the 
Greenfunction are the Fourier transformed causal ones 
from the time ordering. 

We consider the tight-binding Hamiltonian which 
has the matrix structure 



H = 






tj-l 


Vj-i 


tj 



















V J 


tj+l 

















tj+l 


Vj+1 


tj+2 

















tj+2 


Vj+2 


tj+3 



\ 



H ll 


H lr 


H rl 





(A2) 



where we have cut the matrix at the side j . The cut parts 
read separately as 



H rr = (|i>«i(i| + |i+l>*i+i|i> + |i-l>*il<+i; 



i=j+i 

3 



H = 



(\*) v i (*l + K+i>*i+i |») + l*-i)*< |*+i>) 



H lr = \j)t j+1 (j + l\ 
H rl = |j + l)i i+ i(j|. 



(A3) 



Now we employ the general inversion formulas of ma- 
trices composed to 2x2 operators, GB = 1. The upper 
right equation reads 



(A4) 



(A5) 



(A6) 



(A7) 



and from (|A5[) the element G\i follows. The other two 
elements of G are given by interchanging 1 2. 



G11-B12 + G12-B22 — 
from which we obtain 

G12 = — G \\B\iB 22 ■ 
This is used in the upper left equation of GB = 1 

Gii-Bn + G12-B21 = 1 

to find finally 

G11 = (Bn — B12B22 1 B21) 
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These formulas allows to write the Greenfunctions 
(|Aip according to the cutting (|A2[) as 



12 = Lrll-" 



Gn = {z-H ll -H lr G rr H rl ^ 



lr s~ivr 



Gr 1 Tjrl s~ill r < jjrr tjtI fill Tjlr 

21 — ^22^i Lr , O22 — yz — tt —12 Lr ri 



.(A8) 



We consider now the specific energy parts in the diagonal 
Greenfunctions with the help of (|A3[) 



H lr G rr H H = |j) tj + l (J + 1| CT r \ j + 1) tj+l {j\ 



(A9) 



which defines the right-side surface Greenfunction S r 
Analogously we obtain the left-side Greenfunction 



H n G u H ir = \ J + i )tj+1 (j\G u \j)t J+1 (3 + l\ 



(A10) 



These surface Greenfunction obey simple recursion re- 
lations. To see this we cut the left Hamiltonian as 





( - 






H u = 












tj 








-1] 


(!) 








v 3 



) 



(All) 



where H" [j — 1] denotes the left-side Hamiltonian when 
(IA2[) is cutted at j — 1. Due to the equation G n (,z — 
If") = 1 this partition (|A1 1|) enforces the structure of 
the Greenfunction 



G = 



G u [j-1] 



G 



jj 



(A12) 



The left upper part is the Greenfunction G u with the cut 
at j — 1. With the help of the inversion formulas above, 
the right lower element of (I A12[l determines the left-side 
surface Greenfunction 



?l - 



1 



(...0,-*j)G"[7-l] 







y-2 ci 



(A13) 



which establishes the recursion formula for the left-side 
surface Greenfunction. Similiarly one obtains the recur- 
sion for the right-side surface Greenfunction 



S 



1 



Uj+l 



(A14) 



+ 2 



If the tight-binding system is in-between two leads with 
the left material characterized by tj — to and the crystal 



levels Vj = vo for j < and the right material charac- 
terized by tj = ijv+i and Vj — vjv+i for j > N + 1 one 
finds from (|A13j) directly the surface Greenfunction of 
the leads 



s l = — [ z ~ v ° 

to \ 2to 



z - Vp 
2to 



(A15) 



and analogously for Spf + v Together with (|A13|) this de- 
termines all the surface Greenfunctions completely. 

With the knowledge of the surface Greenfunction we 
can now provide the complete Greenfunction (|A7[) . Em- 
ploying the familiar relation A ^ B = i + jj-Bjr + ••■ we 
can write the matrix elements of 

= (k\ G u + G u a \j) (j\ G u + G u a \j) (j\ G u a \j) (j\ G u 
+... |m) 

= Gkm + G kj a C 1 + S j a + ( Sl j a ) 2 + •■■) G ) 
11 , 



-ill 
r jm 



G kj a G l j m 



where we used the abbreviation a 
gously we obtain 



or j.2 



(k\G 12 \m) 
(k\ G 22 \m) 
(k\G 21 \m) 



k jtj+lUj +lr 



{Jr kj+l u '~*j+lr 



i-bs: 

t 2 



(A16) 
Analo- 

(A17) 
(A18) 
(A19) 



with the abbreviation b 

Now we remember the structure of the complete Green- 
function following (|A2|) 



k 
i 

j 

i + i 



V 



■j j + 1- \ 



Gu G12 



G21 G22 



(A20) 



/ 



Comparing G21 of (|A19|) for k = j + 1 and m < j in 
(jAlO) with Gn of (|A16l) for fc = j and m < j in (pO0|) 
we obtain the recursive relation 



Gj+ 



j+l,m 



tj + iSj +1 G^ 



m < j 



(A21) 



which allows to construct the Greenfunction on one side 
from its diagonal part. Analogously we obtain the other 
diagonal part by comparing G22 of (|A18[) for k = j + 1 
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and to > j + 1 in (|A20p with G\i of (|A16|) for k — j and 2. Determine the surface Greenfunctions due to the 

m > j + 1 in (|A20p as recursive relations (|A13jl and (|A14|) 

Gj, m = t j+ xS]G j+hm , m>j + l. (A22) 

The diagonal form is found from Gn of (|A16|1 easily to 1 
be 5i 



|C[l|il ^ Ra? ^ = ,- w -^ - (A25) 

(A23) 



where we used (|A13|) once more. This completes the 
recursive construction of the Greenfunction. 



3. Compute the required elements of the Greenfunc- 
1. Summary of the method tion matrix according to (|A21| ). (|A22[) and (|A23|) 

For computational purposes let us collect the impor- 
tant steps and formulas. 



1. Calculate the surface Greenfunctions of the left and G 



kk 



right lead according to ([Al~5|) z Vk t k+i S k+i t k S k-i 

Gj,m = tj + iS l jGj+i t . m , TO > j + 1 



S r _ 1 ^ z~v N+1 .J^ f z-v N+1 \ \ G j+liin = tj +1 S^ +1 Gj tm , m<j. (A26) 




(A24) 
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